program BGK

  use fft_wrap
  use mpi_wrap
!  use efield
!  use evolve

  use dimensions

  use parameters, only : restart, nt, dt, iDiag, iData, vPhase, runname

  implicit none

  !-- major arrays ---------------------


  real       :: dfu(jv, jz, jy, jx, jDir)       ! distribution function
  real       :: uxMat(jz, jy, jx, jDir)         ! velocity, x-component
  real       :: Ex2D(jz, jy, jx)                ! electric field, x-component 
  real       :: Ez2D(jz, jy, jx)                ! electric field, z-component
  real       :: dens(jz, jy, jx)                ! density

  real       :: filter(jz,jy,jx)                ! conical filter

  real       :: Fbgk(jv,jz)                     ! f_bgk
  real       :: DVbgk(jv,jz)                    ! df_bgk / dv
  real       :: Ebgk(jz)                        ! E_bgk
  real       :: Ubgk(jz)                        ! longitudinal velocity

  complex    :: ftEz(jz, jy, jx)                ! Fourier trasnform of Ez


  !-- variables  ----------------------- 

  real       :: t                               ! time
  integer    :: it                              ! iteration count
  integer    :: iargc                           ! argument count


!-------------------------------------------------------------------------------------

      !--  reading the run_name from the command line

      if(iargc().eq.0) then
         call getarg(0,runname)
         write(*,*)      'Usage: ',trim(runname),' <run name>'
         call mpi_abort(MPI_COMM_WORLD, -1, mpi_err)
         stop
      end if
      call getarg(1,runname)


      !-- set parameters and coordinate space --

      call mpi_wrap_init()

      call read_input()

      call fft_init()

      call set_space_vectors()


      !-- restore from restart file or create initial conditions --

      call maxwell(Fbgk, DVbgk, Ebgk)


      call io_read_bgk(Fbgk)


      !call bgk_test_equil(Fbgk, DVbgk, Ebgk)
      !call mpi_wrap_exit()
      !return

      call efield_bgk(Fbgk, Ebgk)
      call derivV_bgk(Fbgk, DVbgk)
      call sumUz_bgk(Fbgk,  Ubgk)

      Fbgk  = Fbgk / (jDir*npDir)
      DVbgk = DVbgk /(jDir*npDir)

      call set_ic(dfu, uxMat, Fbgk, it)

      t  = it*dt

      call set_filter(filter)

      call evolve_lin_efield( dfu, Ex2D, Ez2D, ftEz, filter)
      call write_bin(dfu, uxMat,  Ex2D, Ez2D, it, it/iData)

      write(msg,*) char(10), ' Begin temporal evolution', char(10)
      call mpi_wrap_msg


!-------------------------------------------------------------------------------------
! BEGIN EVOLUTION
!-------------------------------------------------------------------------------------

  do while (it .le. nt) 

    call evolve_lin_F_long_advect( dfu, dt/2)
    call evolve_lin_F_trans_advect( dfu, Fbgk, uxMat, dt/2)
    call evolve_F_collisions(dfu, dt/2)


    call evolve_lin_efield( dfu, Ex2D, Ez2D, ftEz, filter)
 
    call evolve_lin_F_accel( dfu, Ebgk, DVbgk, Ez2D, dt)

    call evolve_lin_U_long_advect( uxMat, Ubgk, dt)    
    call evolve_lin_U_trans_advect( uxMat, dt)
    call evolve_lin_U_sourceE( uxMat, Ex2D, dt)
 
    call evolve_lin_U_filter( uxMat, filter )


    call evolve_F_collisions(dfu, dt/2)
    call evolve_lin_F_trans_advect( dfu, Fbgk, uxMat, dt/2)
    call evolve_lin_F_long_advect( dfu, dt/2)

    t  = t + dt
    it = it + 1
   

    !-- diagnostics --

    if ( (it/iDiag)*iDiag.eq.it ) then      
   

      call evolve_lin_efield( dfu, Ex2D, Ez2D, ftEz, filter)

      call diagnostics( dfu, Ex2D, Ez2D, uxMat, t)
      write(msg, "(' Diagnostcs, t =' e12.6)") t
      call mpi_wrap_msg

  
      !-- output binary data for restart and post-processing --

      if ( (it/iData)*iData.eq.it ) then

         call write_bin(dfu, uxMat,  Ex2D, Ez2D, it, it/iData)

         write(msg, "(' Output ' i4.4)") it/iData
         call mpi_wrap_msg
       
      end if

    end if

  end do      ! END TEMPORAL LOOP
 

  write(msg,*) char(10), ' End temporal evolution', char(10)
  call mpi_wrap_msg

  !-- write some specrtral spectrum quantities --
       
  call mpi_wrap_exit()
    
  stop;

end


!-------------------------------------------------------------------------------------
subroutine maxwell(Fbgk, DVbgk, Ebgk)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only  :  runname, Vmax, pi, vPhase, iWaveFrame

  implicit none

  real, intent(out)       :: Fbgk(jv,jz)                     ! f_bgk
  real, intent(out)       :: DVbgk(jv,jz)                    ! df_bgk / dv
  real, intent(out)       :: Ebgk(jz)                        ! E_bgk

  real                    :: f0
  integer                 :: iv

  Ebgk(:)=0.0

  ! reminder:  vvec(iv) = -vPhase*iWaveFrame - Vmax + dv*(iv-1)
  ! reminder:  f0vec(:) = ( 1/sqrt(2*pi) )*exp(-0.5*(vvec(:) + vPhase*iWaveFrame)**2)

  do iv = 1,jv
     f0 = ( 1/sqrt(2*pi) )*exp(-0.5*(dv*(iv-1) - Vmax)**2)
     Fbgk(iv,:) = f0
     DVbgk(iv,:) =  -(dv*(iv-1) - Vmax)*f0
  end do

end subroutine maxwell

!-------------------------------------------------------------------------------------
